!-----------------------------------
! A. Gazizov, LNGS, 11.01.11

module coeff
   implicit none

   private AprBD

! ee-pair production on CMB
   integer,parameter,private        &
       :: nee   = 201                 ! number of energy lines	   
   real(8),dimension(nee),private   &
       :: lgee, lgbetaee &	          ! lg(E/eV), lg(beta^{ee}/H_0)
        , dbetaee	  		          ! d\beta^{ee}/dlnE

! N + \gamma -> N' + X on CMB 
   integer,parameter,private        &
       :: necmb = 326                 ! number of energy lines
   real(8),dimension(necmb),private & 
       :: lgecmb                      ! log10(E/eV)

! The case of quick oscillation
! No neutron components are to be taken into account
   real(8), private:: tt   ! foo
   real(8), dimension(necmb), private   &
      :: Ppn                            &
       , Bpp, dBpp, Dpp, dDpp, d2Dpp    & 
       , Bpn, dBpn, Dpn, dDpn, d2Dpn 
   real(8), dimension(13), private:: C

   real(8), dimension(5):: Cour, Cmir	 ! output of this module

!------------------------
         CONTAINS
!------------------------

!-----------------------------------------
   Subroutine ReadPCMB
      implicit none
      integer:: j

! ------------ Arrays for CMB  ---------------
      open(1,file='Data/BBPee.dat',status='old')
         read(1,*)
         read(1,*)
         read(1,'(f8.4,f18.11,E20.11)')               &
            (lgee(j),lgbetaee(j),dbetaee(j), j=1,nee)
      close(1)
!------------- 
     open(1,file="Data/pwij.dat",status='old')
        read(1,'(//(25E15.8))')          & 
       (                                 & 
              lgecmb(j)                  &
            , tt      , Ppn  (j), tt, tt &
            , Bpp  (j), Bpn  (j), tt, tt &
            , Dpp  (j), Dpn  (j), tt, tt &
            , dBpp (j), dBpn (j), tt, tt &
            , dDpp (j), dDpn (j), tt, tt &
            , d2Dpp(j), d2Dpn(j), tt, tt &
            , j=1,necmb                  &   
         )
      close(1)
      print*, ' CMB coefficients are read!'
! ---------- end array for CMB  -------------

   end Subroutine ReadPCMB
!-----------------------------------------

!***************************************************
! Calculates all approximations to the             *
! coefficients  of e^+ + e^- and                   *
! nucleon photoproduction in CMBR  at z=0          *
!***************************************************
!  r = ln(E/eV)                                    *
!***************************************************
! c(1) = Bee, c( 2) = dBee, c( 3) = Ppn                 
! c(4) = Bpp, c( 5) = dBpp, c( 6) = Dpp, c( 7) = dDpp, c( 8) = d2Dpp  
! c(9) = Bpn, c(10) = dBpn, c(11) = Dpn, c(12) = dDpn, c(13) = d2Dpn  

   Subroutine AprBD(vlgE) 
      use msimsl, only: DCSIEZ 
      use commirr, only: r_0, ln10, FLIN 
      implicit none

      integer:: i
      real(8), intent( in):: vlgE
      real(8), dimension(1):: rv     ! lg(E/eV)
      real(8), dimension(13,1):: res
      rv(1) = vlgE


!**********************************************************************
      if(vlgE.le.lgee(1)) then
         res(1,1) = flin(vlgE,lgee(1),lgee(2),lgbetaee(1),lgbetaee(2))  ! lg B^ee
         res(2,1) = flin(vlgE,lgee(1),lgee(2),dbetaee (1),dbetaee (2))  ! dB^ee/dlnE 
      elseif(vlgE.gt.lgee(nee)) then
         res(1,1) = flin(vlgE,lgee(nee-1),lgee(nee),lgbetaee(nee-1),lgbetaee(nee))
         res(2,1) = flin(vlgE,lgee(nee-1),lgee(nee),dbetaee (nee-1),dbetaee (nee))
      else
         CALL DCSIEZ (nee, lgee, lgbetaee,1,rv,res(1,:) )
         CALL DCSIEZ (nee, lgee, dbetaee ,1,rv,res(2,:)	)
      endif
!**********************************************************************
      if(vlgE.le.lgecmb(1)) then
         res( 3,1) = flin(vlgE,lgecmb(1),lgecmb(2),Ppn  (1),Ppn  (2))
         res( 4,1) = flin(vlgE,lgecmb(1),lgecmb(2),Bpp  (1),Bpp  (2))
         res( 5,1) = flin(vlgE,lgecmb(1),lgecmb(2),dBpp (1),dBpp (2))
         res( 6,1) = flin(vlgE,lgecmb(1),lgecmb(2),Dpp  (1),Dpp  (2))
         res( 7,1) = flin(vlgE,lgecmb(1),lgecmb(2),dDpp (1),dDpp (2))
         res( 8,1) = flin(vlgE,lgecmb(1),lgecmb(2),d2Dpp(1),d2Dpp(2))
         res( 9,1) = flin(vlgE,lgecmb(1),lgecmb(2),Bpn  (1),Bpn  (2))
         res(10,1) = flin(vlgE,lgecmb(1),lgecmb(2),dBpn (1),dBpn (2))
         res(11,1) = flin(vlgE,lgecmb(1),lgecmb(2),Dpn  (1),Dpn  (2))
         res(12,1) = flin(vlgE,lgecmb(1),lgecmb(2),dDpn (1),dDpn (2))
         res(13,1) = flin(vlgE,lgecmb(1),lgecmb(2),d2Dpn(1),d2Dpn(2))
      elseif(vlgE.gt.lgecmb(necmb)) then
         res( 3,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),Ppn  (necmb-1),Ppn  (necmb))
         res( 4,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),Bpp  (necmb-1),Bpp  (necmb))
         res( 5,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),dBpp (necmb-1),dBpp (necmb))
         res( 6,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),Dpp  (necmb-1),Dpp  (necmb))
         res( 7,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),dDpp (necmb-1),dDpp (necmb))
         res( 8,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),d2Dpp(necmb-1),d2Dpp(necmb))
         res( 9,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),Bpn  (necmb-1),Bpn  (necmb))
         res(10,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),dBpn (necmb-1),dBpn (necmb))
         res(11,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),Dpn  (necmb-1),Dpn  (necmb))
         res(12,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),dDpn (necmb-1),dDpn (necmb))
         res(13,1) = flin(vlgE,lgecmb(necmb-1),lgecmb(necmb),d2Dpn(necmb-1),d2Dpn(necmb))
      else
         CALL DCSIEZ (necmb, lgecmb,   Ppn ,1, rv, res( 3,:) )
         CALL DCSIEZ (necmb, lgecmb,   Bpp ,1, rv, res( 4,:) )
         CALL DCSIEZ (necmb, lgecmb,  dBpp ,1, rv, res( 5,:) )
         CALL DCSIEZ (necmb, lgecmb,   Dpp ,1, rv, res( 6,:) )
         CALL DCSIEZ (necmb, lgecmb,  dDpp ,1, rv, res( 7,:) )
         CALL DCSIEZ (necmb, lgecmb, d2Dpp ,1, rv, res( 8,:) )
         CALL DCSIEZ (necmb, lgecmb,   Bpn ,1, rv, res( 9,:) )
         CALL DCSIEZ (necmb, lgecmb,  dBpn ,1, rv, res(10,:) )
         CALL DCSIEZ (necmb, lgecmb,   Dpn ,1, rv, res(11,:) )
         CALL DCSIEZ (necmb, lgecmb,  dDpn ,1, rv, res(12,:) )
         CALL DCSIEZ (necmb, lgecmb, d2Dpn ,1, rv, res(13,:) )
      endif

      C(1) = 1d1**res(1,1)
      do i = 2, 13
         C(i) = res(i,1)
      enddo 

   end Subroutine AprBD
!**********************************************************************


!**********************************************************************
!     z1 = 1+z ;      r = ln(E/E_0)
!-----------------------------------------
! Cour(1) - U_r''  I  Cmir(1) - U_m''
! Cour(2) - U_r'   I  Cmir(2) - Um'
! Cour(3) - U_r	   I  Cmir(3) - U_m
! Cour(4) - Ppn_r  I  Cmir(4) - Ppn_m
! Cour(5) - Q_r	   I  Cmir(5) - Q_m
!-----------------------------------------
! c(1) = Bee, c( 2) = dBee, c( 3) = Ppn                 
! c(4) = Bpp, c( 5) = dBpp, c( 6) = Dpp, c( 7) = dDpp, c( 8) = d2Dpp  
! c(9) = Bpn, c(10) = dBpn, c(11) = Dpn, c(12) = dDpn, c(13) = d2Dpn  

   Subroutine AprCoefZ(z1,r)
      use commirr, only: Lambda, Omega_m &
                       , r_0, ln10       &
                       , gam             &
                       , evl             & 
                       , Xmir, Xmir3     &
                       , Osc, r_crit     &
                       , ohne
      implicit none
      integer:: i
      real(8), intent(in):: z1, r
      real(8):: vlgE_our, vlgE_mir, z13, hz, fact_z, RE, REm1, z1evl

      z13    = z1*z1*z1
      z1evl  = z1**evl 
      hz     = dsqrt(Omega_m*z13 + Lambda)
      fact_z = z13/hz

      if(Osc == 'N') then
         RE   = 0d0
         REm1 = 1
      else
         RE = 0.5d0/(1d0 + dexp(2d0*(r - r_crit)))
         REm1 = 1d0 - RE
      endif

!++++++++++++++++++++++  Our World  +++++++++++++++++++++++++++
      vlgE_our  = (r + r_0)/ln10 + dlog10(z1)   ! lg(E/eV) in our world	 

      if( ohne == 'Y' ) then
         Cour(1) = 0d0     
         Cour(2) = 0d0     
         Cour(3) = -2d0     
         Cour(4) = 0d0     
      else      
         call AprBD(vlgE_our)
         do i = 1,13
            C(i) = fact_z*C(i)
         enddo  
         Cour(1) = 0.5d0*( C(6) + REm1*C(11) )
         Cour(2) = 1d0 + C(1) + C(4) + 1.5d0*C( 6) + C( 7)            & 
                 + REm1*(C(9) + 1.5d0*C(11) + C(12) )  
         Cour(3) = -2d0 + C( 1) + C( 2)                               &
                 + C(4) + C( 5) + C( 6) + 1.5d0*C( 7) + 0.5d0*C( 8)   &
                 + REm1*( C(9) + C(10) + C(11) + 1.5d0*C(12)          &
                 + 0.5d0*C(13) )
         Cour(4) = RE*C(3)
      endif
!--- Source generation function ---
! Q = (1+z)^3*(E/E_0)^{-gamma}
!     if(r < 0) then
!         Cour(5) =	fact_z*dexp(-2d0*r)
!     else
         Cour(5) = fact_z*dexp(-gam*r)*z1evl
!     endif
!----------------------------------

!++++++++++++++++  Mirror World  +++++++++++++++
! Here temperature at any epoch z is lower; T(z)_mirr = T(z)_our*Xmir
! The density of mirror photons is lower: n(z)_mirr = n(z)_our*Xmir^3
 
      vlgE_mir  = vlgE_our + dlog10(Xmir)   ! lg(E/eV) in mirror world

      if( ohne == 'Y' ) then
         Cmir(1) = 0d0     
         Cmir(2) = 0d0     
         Cmir(3) = -2d0     
         Cmir(4) = 0d0     
      else      
         call AprBD(vlgE_mir)

         do i = 1,13
            C(i) = fact_z*C(i)*Xmir3
         enddo

         Cmir(1) = 0.5d0*( C(6) + REm1*C(11) )
         Cmir(2) = 1d0 + C(1) + C(4) + 1.5d0*C( 6) + C( 7)            & 
                       + REm1*(C(9) + 1.5d0*C(11) + C(12) )  
         Cmir(3) = -2d0 + C(1) + C(2)                                 &
                 + C(4) + C( 5) + C( 6) + 1.5d0*C( 7) + 0.5d0*C( 8)   &
                 + REm1*( C(9) + C(10) + C(11) + 1.5d0*C(12)          &
                 + 0.5d0*C(13) )
         Cmir(4) = RE*C(3)
      endif
!--- Source generation function ---
! Q = rqqmir*(1+z)^3*(E/E_0)^{-gamma}*(1+z)^evl
!     if(r < 0) then
!         Cmir(5) =	rqqmir*fact_z*dexp(-2d0*r)
!     else
!         Cmir(5) =	rqqmir*fact_z*dexp(-gam*r)*z1evl
!     endif
!----------------------------------
!+++++++++++++++++++++++++++++++++++++++++++++++

   end Subroutine AprCoefZ
!**********************************************************************

end module coeff
